home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMAT6.CPP < prev    next >
C/C++ Source or Header  |  1995-01-17  |  15KB  |  636 lines
  1. //$$ newmat6.cpp            Operators, element access, submatrices
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10.  
  11. //#define REPORT { static ExeCounter ExeCount(__LINE__,6); ++ExeCount; }
  12.  
  13. #define REPORT {}
  14.  
  15. /*************************** general utilities *************************/
  16.  
  17. static int tristore(int n)                      // els in triangular matrix
  18. { return (n*(n+1))/2; }
  19.  
  20.  
  21. /****************************** operators *******************************/
  22.  
  23. Real& Matrix::operator()(int m, int n)
  24. {
  25.    if (m<=0 || m>nrows || n<=0 || n>ncols)
  26.       Throw(IndexException(m,n,*this));
  27.    return store[(m-1)*ncols+n-1];
  28. }
  29.  
  30. Real& SymmetricMatrix::operator()(int m, int n)
  31. {
  32.    if (m<=0 || n<=0 || m>nrows || n>ncols)
  33.       Throw(IndexException(m,n,*this));
  34.    if (m>=n) return store[tristore(m-1)+n-1];
  35.    else return store[tristore(n-1)+m-1];
  36. }
  37.  
  38. Real& UpperTriangularMatrix::operator()(int m, int n)
  39. {
  40.    if (m<=0 || n<m || n>ncols)
  41.       Throw(IndexException(m,n,*this));
  42.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  43. }
  44.  
  45. Real& LowerTriangularMatrix::operator()(int m, int n)
  46. {
  47.    if (n<=0 || m<n || m>nrows)
  48.       Throw(IndexException(m,n,*this));
  49.    return store[tristore(m-1)+n-1];
  50. }
  51.  
  52. Real& DiagonalMatrix::operator()(int m, int n)
  53. {
  54.    if (n<=0 || m!=n || m>nrows || n>ncols)
  55.       Throw(IndexException(m,n,*this));
  56.    return store[n-1];
  57. }
  58.  
  59. Real& DiagonalMatrix::operator()(int m)
  60. {
  61.    if (m<=0 || m>nrows) Throw(IndexException(m,*this));
  62.    return store[m-1];
  63. }
  64.  
  65. Real& ColumnVector::operator()(int m)
  66. {
  67.    if (m<=0 || m> nrows) Throw(IndexException(m,*this));
  68.    return store[m-1];
  69. }
  70.  
  71. Real& RowVector::operator()(int n)
  72. {
  73.    if (n<=0 || n> ncols) Throw(IndexException(n,*this));
  74.    return store[n-1];
  75. }
  76.  
  77. Real& BandMatrix::operator()(int m, int n)
  78. {
  79.    int w = upper+lower+1; int i = lower+n-m;
  80.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  81.       Throw(IndexException(m,n,*this));
  82.    return store[w*(m-1)+i];
  83. }
  84.  
  85. Real& UpperBandMatrix::operator()(int m, int n)
  86. {
  87.    int w = upper+1; int i = n-m;
  88.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  89.       Throw(IndexException(m,n,*this));
  90.    return store[w*(m-1)+i];
  91. }
  92.  
  93. Real& LowerBandMatrix::operator()(int m, int n)
  94. {
  95.    int w = lower+1; int i = lower+n-m;
  96.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  97.       Throw(IndexException(m,n,*this));
  98.    return store[w*(m-1)+i];
  99. }
  100.  
  101. Real& SymmetricBandMatrix::operator()(int m, int n)
  102. {
  103.    int w = lower+1;
  104.    if (m>=n)
  105.    {
  106.       int i = lower+n-m;
  107.       if ( m>nrows || n<=0 || i<0 )
  108.          Throw(IndexException(m,n,*this));
  109.       return store[w*(m-1)+i];
  110.    }
  111.    else
  112.    {
  113.       int i = lower+m-n;
  114.       if ( n>nrows || m<=0 || i<0 )
  115.          Throw(IndexException(m,n,*this));
  116.       return store[w*(n-1)+i];
  117.    }
  118. }
  119.  
  120. #ifndef __ZTC__
  121.  
  122. Real Matrix::operator()(int m, int n) const
  123. {
  124.    if (m<=0 || m>nrows || n<=0 || n>ncols)
  125.       Throw(IndexException(m,n,*this));
  126.    return store[(m-1)*ncols+n-1];
  127. }
  128.  
  129. Real SymmetricMatrix::operator()(int m, int n) const
  130. {
  131.    if (m<=0 || n<=0 || m>nrows || n>ncols)
  132.       Throw(IndexException(m,n,*this));
  133.    if (m>=n) return store[tristore(m-1)+n-1];
  134.    else return store[tristore(n-1)+m-1];
  135. }
  136.  
  137. Real UpperTriangularMatrix::operator()(int m, int n) const
  138. {
  139.    if (m<=0 || n<m || n>ncols)
  140.       Throw(IndexException(m,n,*this));
  141.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  142. }
  143.  
  144. Real LowerTriangularMatrix::operator()(int m, int n) const
  145. {
  146.    if (n<=0 || m<n || m>nrows)
  147.       Throw(IndexException(m,n,*this));
  148.    return store[tristore(m-1)+n-1];
  149. }
  150.  
  151. Real DiagonalMatrix::operator()(int m, int n) const
  152. {
  153.    if (n<=0 || m!=n || m>nrows || n>ncols)
  154.       Throw(IndexException(m,n,*this));
  155.    return store[n-1];
  156. }
  157.  
  158. Real DiagonalMatrix::operator()(int m) const
  159. {
  160.    if (m<=0 || m>nrows) Throw(IndexException(m,*this));
  161.    return store[m-1];
  162. }
  163.  
  164. Real ColumnVector::operator()(int m) const
  165. {
  166.    if (m<=0 || m> nrows) Throw(IndexException(m,*this));
  167.    return store[m-1];
  168. }
  169.  
  170. Real RowVector::operator()(int n) const
  171. {
  172.    if (n<=0 || n> ncols) Throw(IndexException(n,*this));
  173.    return store[n-1];
  174. }
  175.  
  176. Real BandMatrix::operator()(int m, int n) const
  177. {
  178.    int w = upper+lower+1; int i = lower+n-m;
  179.    if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  180.       Throw(IndexException(m,n,*this));
  181.    return store[w*(m-1)+i];
  182. }
  183.  
  184. Real SymmetricBandMatrix::operator()(int m, int n) const
  185. {
  186.    int w = lower+1;
  187.    if (m>=n)
  188.    {
  189.       int i = lower+n-m;
  190.       if ( m>nrows || n<=0 || i<0 )
  191.          Throw(IndexException(m,n,*this));
  192.       return store[w*(m-1)+i];
  193.    }
  194.    else
  195.    {
  196.       int i = lower+m-n;
  197.       if ( n>nrows || m<=0 || i<0 )
  198.          Throw(IndexException(m,n,*this));
  199.       return store[w*(n-1)+i];
  200.    }
  201. }
  202.  
  203. #endif
  204.  
  205. Real BaseMatrix::AsScalar() const
  206. {
  207.    REPORT
  208.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  209.    if (gm->nrows!=1 || gm->ncols!=1)
  210.    {
  211.       Try
  212.       {
  213.          Tracer tr("AsScalar");
  214.      Throw(ProgramException("Cannot convert to scalar", *gm));
  215.       }
  216.       CatchAll { gm->tDelete(); Bounce; }
  217.    }
  218.    Real x = *(gm->store); gm->tDelete(); return x;
  219. }
  220.  
  221. #ifdef TEMPS_DESTROYED_QUICKLY
  222.  
  223. AddedMatrix& BaseMatrix::operator+(const BaseMatrix& bm) const
  224. {
  225.    REPORT
  226.    AddedMatrix* x = new AddedMatrix(this, &bm);
  227.    MatrixErrorNoSpace(x);
  228.    return *x;
  229. }
  230.  
  231. SPMatrix& SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
  232. {
  233.    REPORT
  234.    SPMatrix* x = new SPMatrix(&bm1, &bm2);
  235.    MatrixErrorNoSpace(x);
  236.    return *x;
  237. }
  238.  
  239. MultipliedMatrix& BaseMatrix::operator*(const BaseMatrix& bm) const
  240. {
  241.    REPORT
  242.    MultipliedMatrix* x = new MultipliedMatrix(this, &bm);
  243.    MatrixErrorNoSpace(x);
  244.    return *x;
  245. }
  246.  
  247. ConcatenatedMatrix& BaseMatrix::operator|(const BaseMatrix& bm) const
  248. {
  249.    REPORT
  250.    ConcatenatedMatrix* x = new ConcatenatedMatrix(this, &bm);
  251.    MatrixErrorNoSpace(x);
  252.    return *x;
  253. }
  254.  
  255. StackedMatrix& BaseMatrix::operator&(const BaseMatrix& bm) const
  256. {
  257.    REPORT
  258.    StackedMatrix* x = new StackedMatrix(this, &bm);
  259.    MatrixErrorNoSpace(x);
  260.    return *x;
  261. }
  262.  
  263. //SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) const
  264. SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx)
  265. {
  266.    REPORT
  267.    SolvedMatrix* x;
  268.    Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
  269.    CatchAll { delete this; Bounce; }
  270.    delete this;                // since we are using bm rather than this
  271.    return *x;
  272. }
  273.  
  274. SubtractedMatrix& BaseMatrix::operator-(const BaseMatrix& bm) const
  275. {
  276.    REPORT
  277.    SubtractedMatrix* x = new SubtractedMatrix(this, &bm);
  278.    MatrixErrorNoSpace(x);
  279.    return *x;
  280. }
  281.  
  282. ShiftedMatrix& BaseMatrix::operator+(Real f) const
  283. {
  284.    REPORT
  285.    ShiftedMatrix* x = new ShiftedMatrix(this, f);
  286.    MatrixErrorNoSpace(x);
  287.    return *x;
  288. }
  289.  
  290. ScaledMatrix& BaseMatrix::operator*(Real f) const
  291. {
  292.    REPORT
  293.    ScaledMatrix* x = new ScaledMatrix(this, f);
  294.    MatrixErrorNoSpace(x);
  295.    return *x;
  296. }
  297.  
  298. ScaledMatrix& BaseMatrix::operator/(Real f) const
  299. {
  300.    REPORT
  301.    ScaledMatrix* x = new ScaledMatrix(this, 1.0/f);
  302.    MatrixErrorNoSpace(x);
  303.    return *x;
  304. }
  305.  
  306. ShiftedMatrix& BaseMatrix::operator-(Real f) const
  307. {
  308.    REPORT
  309.    ShiftedMatrix* x = new ShiftedMatrix(this, -f);
  310.    MatrixErrorNoSpace(x);
  311.    return *x;
  312. }
  313.  
  314. TransposedMatrix& BaseMatrix::t() const
  315. {
  316.    REPORT
  317.    TransposedMatrix* x = new TransposedMatrix(this);
  318.    MatrixErrorNoSpace(x);
  319.    return *x;
  320. }
  321.  
  322. NegatedMatrix& BaseMatrix::operator-() const
  323. {
  324.    REPORT
  325.    NegatedMatrix* x = new NegatedMatrix(this);
  326.    MatrixErrorNoSpace(x);
  327.    return *x;
  328. }
  329.  
  330. InvertedMatrix& BaseMatrix::i() const
  331. {
  332.    REPORT
  333.    InvertedMatrix* x = new InvertedMatrix(this);
  334.    MatrixErrorNoSpace(x);
  335.    return *x;
  336. }
  337.  
  338. ConstMatrix& GeneralMatrix::c() const
  339. {
  340.    if (tag != -1)
  341.       Throw(ProgramException(".c() applied to temporary matrix"));
  342.    REPORT
  343.    ConstMatrix* x = new ConstMatrix(this);
  344.    MatrixErrorNoSpace(x);
  345.    return *x;
  346. }
  347.  
  348. RowedMatrix& BaseMatrix::AsRow() const
  349. {
  350.    REPORT
  351.    RowedMatrix* x = new RowedMatrix(this);
  352.    MatrixErrorNoSpace(x);
  353.    return *x;
  354. }
  355.  
  356. ColedMatrix& BaseMatrix::AsColumn() const
  357. {
  358.    REPORT
  359.    ColedMatrix* x = new ColedMatrix(this);
  360.    MatrixErrorNoSpace(x);
  361.    return *x;
  362. }
  363.  
  364. DiagedMatrix& BaseMatrix::AsDiagonal() const
  365. {
  366.    REPORT
  367.    DiagedMatrix* x = new DiagedMatrix(this);
  368.    MatrixErrorNoSpace(x);
  369.    return *x;
  370. }
  371.  
  372. MatedMatrix& BaseMatrix::AsMatrix(int nrx, int ncx) const
  373. {
  374.    REPORT
  375.    MatedMatrix* x = new MatedMatrix(this,nrx,ncx);
  376.    MatrixErrorNoSpace(x);
  377.    return *x;
  378. }
  379.  
  380. #else
  381.  
  382. AddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
  383. { REPORT return AddedMatrix(this, &bm); }
  384.  
  385. SPMatrix SP(const BaseMatrix& bm1,const BaseMatrix& bm2)
  386. { REPORT return SPMatrix(&bm1, &bm2); }
  387.  
  388. MultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
  389. { REPORT return MultipliedMatrix(this, &bm); }
  390.  
  391. ConcatenatedMatrix BaseMatrix::operator|(const BaseMatrix& bm) const
  392. { REPORT return ConcatenatedMatrix(this, &bm); }
  393.  
  394. StackedMatrix BaseMatrix::operator&(const BaseMatrix& bm) const
  395. { REPORT return StackedMatrix(this, &bm); }
  396.  
  397. SolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
  398. { REPORT return SolvedMatrix(bm, &bmx); }
  399.  
  400. SubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
  401. { REPORT return SubtractedMatrix(this, &bm); }
  402.  
  403. ShiftedMatrix BaseMatrix::operator+(Real f) const
  404. { REPORT return ShiftedMatrix(this, f); }
  405.  
  406. ScaledMatrix BaseMatrix::operator*(Real f) const
  407. { REPORT return ScaledMatrix(this, f); }
  408.  
  409. ScaledMatrix BaseMatrix::operator/(Real f) const
  410. { REPORT return ScaledMatrix(this, 1.0/f); }
  411.  
  412. ShiftedMatrix BaseMatrix::operator-(Real f) const
  413. { REPORT return ShiftedMatrix(this, -f); }
  414.  
  415. TransposedMatrix BaseMatrix::t() const
  416. { REPORT return TransposedMatrix(this); }
  417.  
  418. NegatedMatrix BaseMatrix::operator-() const
  419. { REPORT return NegatedMatrix(this); }
  420.  
  421. InvertedMatrix BaseMatrix::i() const
  422. { REPORT return InvertedMatrix(this); }
  423.  
  424. ConstMatrix GeneralMatrix::c() const
  425. {
  426.    if (tag != -1)
  427.       Throw(ProgramException(".c() applied to temporary matrix"));
  428.    REPORT return ConstMatrix(this);
  429. }
  430.  
  431. RowedMatrix BaseMatrix::AsRow() const
  432. { REPORT return RowedMatrix(this); }
  433.  
  434. ColedMatrix BaseMatrix::AsColumn() const
  435. { REPORT return ColedMatrix(this); }
  436.  
  437. DiagedMatrix BaseMatrix::AsDiagonal() const
  438. { REPORT return DiagedMatrix(this); }
  439.  
  440. MatedMatrix BaseMatrix::AsMatrix(int nrx, int ncx) const
  441. { REPORT return MatedMatrix(this,nrx,ncx); }
  442.  
  443. #endif
  444.  
  445. void GeneralMatrix::operator=(Real f)
  446. { REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
  447.  
  448. void Matrix::operator=(const BaseMatrix& X)
  449. {
  450.    REPORT //CheckConversion(X);
  451.    MatrixConversionCheck mcc;
  452.    Eq(X,MatrixType::Rt);
  453. } 
  454.  
  455. void RowVector::operator=(const BaseMatrix& X)
  456. {
  457.    REPORT  // CheckConversion(X);
  458.    MatrixConversionCheck mcc;
  459.    Eq(X,MatrixType::RV);
  460.    if (nrows!=1)
  461.    {
  462.       Tracer tr("RowVector(=)");
  463.       Throw(VectorException(*this));
  464.    }
  465. }
  466.  
  467. void ColumnVector::operator=(const BaseMatrix& X)
  468. {
  469.    REPORT //CheckConversion(X);
  470.    MatrixConversionCheck mcc;
  471.    Eq(X,MatrixType::CV);
  472.    if (ncols!=1)
  473.    {
  474.       Tracer tr("ColumnVector(=)");
  475.       Throw(VectorException(*this));
  476.    }
  477. }
  478.  
  479. void SymmetricMatrix::operator=(const BaseMatrix& X)
  480. {
  481.    REPORT // CheckConversion(X);
  482.    MatrixConversionCheck mcc;
  483.    Eq(X,MatrixType::Sm);
  484. }
  485.  
  486. void UpperTriangularMatrix::operator=(const BaseMatrix& X)
  487. {
  488.    REPORT //CheckConversion(X);
  489.    MatrixConversionCheck mcc;
  490.    Eq(X,MatrixType::UT);
  491. }
  492.  
  493. void LowerTriangularMatrix::operator=(const BaseMatrix& X)
  494. {
  495.    REPORT //CheckConversion(X);
  496.    MatrixConversionCheck mcc;
  497.    Eq(X,MatrixType::LT);
  498. }
  499.  
  500. void DiagonalMatrix::operator=(const BaseMatrix& X)
  501. {
  502.    REPORT // CheckConversion(X);
  503.    MatrixConversionCheck mcc;
  504.    Eq(X,MatrixType::Dg);
  505. }
  506.  
  507. void GeneralMatrix::operator<<(const Real* r)
  508. {
  509.    REPORT
  510.    int i = storage; Real* s=store;
  511.    while(i--) *s++ = *r++;
  512. }
  513.  
  514.  
  515. void GenericMatrix::operator=(const GenericMatrix& bmx)
  516. {
  517.    if (&bmx != this) { REPORT if (gm) delete gm; gm = bmx.gm->Image();}
  518.    else { REPORT }
  519.    gm->Protect(); 
  520. }
  521.  
  522. void GenericMatrix::operator=(const BaseMatrix& bmx)
  523. {
  524.    if (gm)
  525.    {
  526.       int counter=bmx.search(gm);
  527.       if (counter==0) { REPORT delete gm; gm=0; }
  528.       else { REPORT gm->Release(counter); }
  529.    }
  530.    else { REPORT }
  531.    GeneralMatrix* gmx = ((BaseMatrix&)bmx).Evaluate();
  532.    if (gmx != gm) { REPORT if (gm) delete gm; gm = gmx->Image(); }
  533.    else { REPORT }
  534.    gm->Protect(); 
  535. }
  536.  
  537. /************************* element access *********************************/
  538.  
  539. Real& Matrix::element(int m, int n)
  540. {
  541.    if (m<0 || m>= nrows || n<0 || n>= ncols)
  542.       Throw(IndexException(m,n,*this,TRUE));
  543.    return store[m*ncols+n];
  544. }
  545.  
  546. Real& SymmetricMatrix::element(int m, int n)
  547. {
  548.    if (m<0 || n<0 || m >= nrows || n>=ncols)
  549.       Throw(IndexException(m,n,*this,TRUE));
  550.    if (m>=n) return store[tristore(m)+n];
  551.    else return store[tristore(n)+m];
  552. }
  553.  
  554. Real& UpperTriangularMatrix::element(int m, int n)
  555. {
  556.    if (m<0 || n<m || n>=ncols)
  557.       Throw(IndexException(m,n,*this,TRUE));
  558.    return store[m*ncols+n-tristore(m)];
  559. }
  560.  
  561. Real& LowerTriangularMatrix::element(int m, int n)
  562. {
  563.    if (n<0 || m<n || m>=nrows)
  564.       Throw(IndexException(m,n,*this,TRUE));
  565.    return store[tristore(m)+n];
  566. }
  567.  
  568. Real& DiagonalMatrix::element(int m, int n)
  569. {
  570.    if (n<0 || m!=n || m>=nrows || n>=ncols)
  571.       Throw(IndexException(m,n,*this,TRUE));
  572.    return store[n];
  573. }
  574.  
  575. Real& DiagonalMatrix::element(int m)
  576. {
  577.    if (m<0 || m>=nrows) Throw(IndexException(m,*this,TRUE));
  578.    return store[m];
  579. }
  580.  
  581. Real& ColumnVector::element(int m)
  582. {
  583.    if (m<0 || m>= nrows) Throw(IndexException(m,*this,TRUE));
  584.    return store[m];
  585. }
  586.  
  587. Real& RowVector::element(int n)
  588. {
  589.    if (n<0 || n>= ncols)  Throw(IndexException(n,*this,TRUE));
  590.    return store[n];
  591. }
  592.  
  593. Real& BandMatrix::element(int m, int n)
  594. {
  595.    int w = upper+lower+1; int i = lower+n-m;
  596.    if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  597.       Throw(IndexException(m,n,*this,TRUE));
  598.    return store[w*m+i];
  599. }
  600.  
  601. Real& UpperBandMatrix::element(int m, int n)
  602. {
  603.    int w = upper+1; int i = n-m;
  604.    if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  605.       Throw(IndexException(m,n,*this,TRUE));
  606.    return store[w*m+i];
  607. }
  608.  
  609. Real& LowerBandMatrix::element(int m, int n)
  610. {
  611.    int w = lower+1; int i = lower+n-m;
  612.    if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  613.       Throw(IndexException(m,n,*this,TRUE));
  614.    return store[w*m+i];
  615. }
  616.  
  617. Real& SymmetricBandMatrix::element(int m, int n)
  618. {
  619.    int w = lower+1;
  620.    if (m>=n)
  621.    {
  622.       int i = lower+n-m;
  623.       if ( m>=nrows || n<0 || i<0 )
  624.          Throw(IndexException(m,n,*this,TRUE));
  625.       return store[w*m+i];
  626.    }
  627.    else
  628.    {
  629.       int i = lower+m-n;
  630.       if ( n>=nrows || m<0 || i<0 )
  631.          Throw(IndexException(m,n,*this,TRUE));
  632.       return store[w*n+i];
  633.    }
  634. }
  635.  
  636.